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Abstract 

We use numerical simulations to probe the dynamics of concentrated suspensions of spherical 
microswimmers interacting hydrodynamically. Previous work in the dilute limit predicted orienta- 
tional instabilities of aligned suspensions for both pusher and puller swimmers, which we confirm 
computationally. Unlike previous work, we show that isotropic suspensions of spherical swimmers 
are also always unstable. Both types of initial conditions develop long-time polar order, of a na- 
ture which depends on the hydrodynamic signature of the swimmer but very weakly on the volume 
fraction up to very high volume fractions. 



Living fluids and chemically active colloidal dispersions are modern examples of nonequi- 
librium hydrodynamic phenomena, presenting tantalizing avenues for both research and 
industrial application. Recent experiments on motile particles, from collections of microor- 
ganisms m El El [211 ESI [30] to self-propelled colloids [HI [9l [IH [22j , exhibit pattern forming 
behavior and enhanced transport characteristics that pose fundamental questions for the 
nonequilibrium statistical mechanics of active systems [20j, and have implications in bio- 
and nano-engineering [29j. One area of particular interest concerns the hydrodynamics of 
microorganisms swimming in viscous fluids at zero Reynolds number [3l [131 [E] • 

Several methodologies have been developed in the past to address the emergence of col- 
lective locomotion, corresponding to either microscopic or macroscopic formulations. Active 
hydrodynamic equations developed from non-equilibrium kinetic theory has been the pre- 
vailing microscopic approach to the system [U [23l [26]. By modeling microswimmers as 
force-dipoles, these theories build continuum dynamical equations for flelds quantifying the 
long- wavelength properties of suspensions of self-propelled particles. The difliculty in deal- 
ing with interacting particles leads to the necessity of a dilute assumption, and the lack of 
any specifled structure to the microswimmers (beyond oriented point singularities). 

At the other end of the spectrum, thermodynamic models of active media give access 
to nonlinearities and coupled modes that cannot be derived from a dilute formulation due 
to relevant (allowed) terms in a dynamical equation [I6l[28]. These "flocking" models lead 
to rich dynamics that cannot be captured by the microscopic models, at the expense of 
introducing phenomenological parameters that may be diflicult to explain physically and 
derive from microscopic considerations [20j. 

Here we use a model spherical microswimmer to address orientational order computa- 
tionally. With our approach both semi-dilute and concentrated suspensions of spherical 
swimmers can be considered. Many microorganisms, such as B. subtilis or E. coli are elon- 
gated swimmers, and as a result much work has been done to model the locomotion of rodlike 
active particles, as well as fabricate synthetic devices of similar aspect ratio; however, col- 
loidal microswimmers, active oil droplets or self-propelled vesicles are likely to be spherical 
in shape, not to mention the spherical organisms like Volvox which appear readily in nature. 
For all of these reasons it is important to develop an understanding of the hydrodynamic in- 
teractions between active spherical particles. Contrary to predictions from dilute continuum 
theories, we flnd that isotropic suspensions of spherical swimmers are unstable, and evolve 
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FIG. 1: Flow streamlines of isolated squirmers in the swimming frame (top; a to c) and lab frame 
(bottom; d to f). Left (a and d): Pusher with a negative stresslet (/3 = —5). Center (b and e): 
Potential flow developed by a squirmer with /3 = 0. Right (c and f): Puller with a positive stresslet 
(/3 = +5). 

dynamically to take on a long-time state of polar order. This state of polar order, which 
exists in phenomenological flocking models, is thus shown here to arise from hydrodynamic 
interactions. 

We simulate a system of spherical swimmers in a cubic box of volume L^, where L is 
determined from the number of swimmers and preset volume fraction 0, i.e. (j) = 4:7rN/{3L^) 
(the radius of the swimmer is 1) and periodic boundary conditions. The swimmer we use, 
termed a squirmer [21 [T5j, is a spherical particle that has a prescribed axisymmetric tan- 
gential velocity distribution on its surface. We impose uo{9) = BiPi{co89) + B2P2{cos9)^ 
where Pn is the n^^ Legendre polynomial, and 9 = deflnes the direction e in which the 
squirmer swims, with speed (for a solitary squirmer) U ^ Bi {U is set to 1). Fluid dis- 
turbances in the far fleld are governed by the "stresslet" (or force-dipole) of the organisms, 
quantifled by the dimensionless quantity /3 = B2/B1. The flow fleld decays as ^ 
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FIG. 2: A snapshot of a simulation with (/) = 0.1, /3 = 1 (pullers), and = 64. The computational 
cell is in the middle, with identical copies surrounding this periodic box. Dark dots on the squirmers 
represent the rear (i.e. = tt), while light dots represent the swimming direction, = 0. 

far from the swimmer, and, in the absence of thermal fluctuations, the stresslet dominates 
the long-range interactions p!3], as illustrated in Fig. [l] Some microorganisms, like the al- 
gae Chlamydomonas^ generate thrust in front of their bodies pulling themselves therefore 
through the fluid ( "pullers" , with /5 > 0) while most flagellated cells such the bacteria E. coli^ 
or spermatozoa, generate thrust behind them, and instead are being pushed from the back 
("pushers", with /3 < 0). Typical values for /3 range from ?^ — 1 for bacteria like E. Coli [6j, 
^ for Volvox and artiflcially created squirmers |l7l EU [27] , to ^ +1 for Chlamydomonas 

The swimming kinematics of the = 64 squirmers are calculated by enforcing instanta- 
neously the condition of force- and torque-free swimming. The computational approach is 
based on Stokesian dynamics, with an analytical treatment of lubrication forces for closely 
separated swimmers, as well as short-range repulsive forces to prevent particle overlap, as de- 
scribed in detail in Refs. JIl[l2]. A typical snapshot of the simulation is displayed in Fig.[2| 
The two parameters characterizing the collective swimming dynamics are thus the swimmer 
volume fraction, 0, and the swimmer stresslet, /3. In order to probe the development of 
order in our system, we deflne a polar order parameter, P, based on the orientation vector 
e of the particles, namely P{t) = \ J2f If every particle is swimming in the same 

direction (polar order) then P = 1, while for isotropic orientation we expect P ^ l/^/N. 

In the dilute limit, continuum theories for slender self-propelled rods have predicted that 
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FIG. 3: Global order of semi-dilute suspensions of pushers and pullers starting from aligned 
state, (a): Order parameter, for pushers (/3 = —1), pullers (/3 = 1), and potential swimmers 

(/3 = 0). Each simulation starts from an initially aligned state (with random positions) and decays 
to a finite order, Poo, after a characteristic decay time that depends on the stresslet coefficient /3. 
(b): Long-time order parameter. Poo? foi" initially aligned suspensions. Each data point represents 
an average over an ensemble (five separate realizations) and time (performed after initial transient). 
Error bars represent standard deviations of the ensemble averaging process. The horizontal dashed 
line indicates the expected result for isotropic suspensions [N — 64). 

aligned suspensions are always unstable, for both pushers and pullers. On the contrary, 
isotropically oriented states are only unstable for pushers (/5 < 0) while pullers remain 
in an isotropic state [23j. The physical nature of the isotropic instability is described as 
resulting from long-range hydrodynamic extensional disturbances that cause reorientation 
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of anisotropic particles, for which spherical swimmers are immune. In the dilute limit, 
no instability of isotropic suspensions due to long-range hydrodynamic interactions is thus 
expected to occur for spherical particles [23l [26] . 

In Fig. [3^ we plot the time-evolution of the polar order parameter, for suspensions 

starting in the aligned state with three different stresslet values, /3 = —1 (pusher), (po- 
tential swimmer) and +1 (puller). The initial positions of the swimmers are taken to be 
random. Over the semi-dilute to concentrated range of volume fractions, = 0.1 — 0.5 (the 
results in Fig. [3^ are shown for = 0.1), we observe the system to systematically decay from 
perfect order (P = 1) to some finite long-time value (0 < Poo < !)• The decay time over 
which the suspension is driven to this new ordered state depends on the stresslet value /5, 
and larger values of the stresslet disturb the fiuid more violently, causing reorientation, and 
loss of polar order, more quickly; /3 can thus be interpreted as the speed at which orientation 
decorrelation propagates throughout the suspension. 

To further characterize long-time global order, we perform ensemble averages on five 
realizations of our simulations, starting with an aligned orientations and random positions, 
for volume fractions in the range = 0.1 to 0.5 and stresslet s varying from /3 = —2 to 
+2. Upon doing ensemble averages, we define a long-time order parameter Poo by averaging 
over the time period after the initial decay from alignment until the end of the simulation; 
the results are shown in Fig. [8)3. In the range of volume fractions studied, the long-time 
order is essentially independent of the value of 0, but strongly depends on the value of 
the stresslet /3. Larger values of lead to increased swimmer-swimmer reorientations due 
hydrodynamic interactions, and thus lead to decreased values of Poo- In addition, we observe 
a marked asymmetry between pushers and pullers, and for a given value of pullers are 
systematically more ordered. For large values of the global order disappears, as shown 
by the dashed line in Fig. denoting the isotropy that is expected in our system. 

This is in contrast with analytical predictions in the dilute limit where suspensions of pullers 
are expected to be systematically driven to orientational isotropy [23] . 

It has been suggested in the past that near-field interactions are the dominant mechanism 
for polar ordering [llj. To address the pusher-puller asymmetry, we take a detailed look at 
the fiow field near the swimmers. The far-field difference between a pusher and a puller leads 
in the near-field to a change in the location of the stagnation point, see Fig. [T] for a pusher 
this stagnation point leads the swimmer, while for a puller it trails it. This asymmetry causes 
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FIG. 4: Time-evolution of the order parameter, P{t), for initially aligned (top) and isotropic 
(bottom) suspensions (in both cases initial positions are random). The case shown is for a puller 
swimmer (/3 = 1) with a volume fraction (j) = 0.01. Both aligned and isotropic states are unstable 
and approach similar finale state of global orient ational order at long times. 

the different coUisional situations (head to head, head to tail, and tail to tail) to produce 
different reorientations between pushers and pullers. From previous numerical work, it is 
known that head-to-head interactions are far more likely to occur [lOj. For pushers the head- 
to-head orientation is stable to rotation (see Fig. [l]l): the presence of a stagnation point 
establishes vortices near the surface of the swimmer, and in a head to head collision both 
pairs of vortices interact in a manner as to maintain this orientation. Pullers, on the other 
hand, are unstable in this configuration, as their vortices trail in such a collision. In terms 
of contributing to order, head to head orientations yields P ^ 0, and thus any instability of 
this configuration, as expected for pullers, will lead to an increase of polar order, hence the 
asymmetry between pushers and pullers seen in Fig. [3)3. 

In contrast to the aligned case, orientation instabilities in isotropic suspensions was pre- 
dicted by continuum theories to exist only for pushers, and only then when the swimmers 
have a nonspherical shape [23j . We performed numerical simulations to probe the long-time 
behavior of initially isotropic suspensions. We show in Fig. |4]the evolution in time of the 
order parameter for suspensions of pullers (/3 = 1) at very low volume fraction (0 = 0.01), 
for both aligned (top) and isotropic (bottom) initial conditions. The results for pushers are 
similar. We see not only that both sets of initial conditions are unstable, but also that both 
are driven at long times to a similar, and intermediate, value of the order parameter. The 
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FIG. 5: Polar order as function of volume fraction, (a) Samples of order decay for initially-aligned 
pullers (/3 = 1); (b) Long-time order as a function of volume fraction, for pullers (/3 = 1). For 
semi-dilute suspensions there is only a weak dependence on the volume fraction. At larger volume 
fractions the order drops sharply to isotropy (indicated by dashed line). Error bars show standard 
deviation over many realizations. 

time scale over which the instability takes place is long compared to the case of more con- 
centrated suspensions, and thus it is plausible that in the limit ^ ^ the instability would 
disappear, allowing us to reconcile our numerical results with those of continuum theories. 

How does increasing the volume fraction affect this order? We show in Fig. [5^ the 
evolution in time of the order parameter for different volume fractions, 0, of pullers (/3 = 1). 
The averaged value of the order parameter. Poo, is plotted as a function of the volume fraction 
in Fig. ^jp. We see that as the volume fraction starts increasing away from the dilute limit. 
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FIG. 6: Behavior at high volume fraction for /3 = 0. (a): Dynamics starting from two different 
initial positions for = 0.5; the fluctuations about the long time average Pqo become large at 
high volume fractions; (b): Dynamics starting from three different initial positions for (f) = 0.6; 
above ^ 0.5, the transition to isotropy can show significant delay times, likely associated with 
an inability to sample all of phase space. 

the average value of the order parameter changes only slightly over a wide range of volume 
fractions. When the volume fraction reaches ^ 0.5, fluctuations in the order parameter 
are observed to become larger, until finally the order disappears near ^ 0.6. Suspensions 
of pushers display a similar dependence, although with a decreased overall magnitude of 
the order parameter due to the pusher-puller asymmetry noted earlier. The fluctuating 
behavior observed at high concentrations is further illustrated in Fig. [6^. We also show in 
Fig. (6)3 how different initial conditions for the same dense suspension can yield drastically 
different results, indicating that particles are "trapped" in their initial conditions before 
finally escaping to the isotropic steady state. 

In summary, we have used in this paper numerical simulations to address the instabilities 
and long-time order of semi-dilute and dense suspensions of spherical swimmers. We have 
shown that spherical squirmers, starting from either an aligned or isotropic state, develop 
long-time polar order due to hydrodynamic interactions, of a nature which depends on the 
hydrodynamic signature of the swimmer but very weakly on the volume fraction up to the 
dense regime. Our results show thus non-trivial differences with dilute, dipole continuum 
models, but display similarities with phenomenological flocking models. 

This work was funded in part by the National Science Foundation through grant CBET- 
0746285 (EL) and a 2010 East Asia and Pacific Summer Institute Research Fellowship 
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